#!/usr/bin/env python
# -*- coding: utf-8 -*-

""" By Martin Senande-Rivera
	For Spatial and temporal expansion of global wildland fire activity in response to climate change """

import numpy as np
import pandas as pd
from pandas import Series, DataFrame
import xarray as xr
import datetime as dt
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 14})
plt.rcParams.update({'lines.linewidth': 2})
from matplotlib.ticker import (MultipleLocator, 
                               FormatStrFormatter, 
                               AutoMinorLocator) 

### File reading ###
####################

ruta_gfed4='../../DATA/GFED4/'	# GFED4 data path
ds_ba_a = xr.open_dataset(ruta_gfed4+'GFED4_BA_a_mean.nc')
ds_ba_m = xr.open_dataset(ruta_gfed4+'GFED4_BA_m_mean.nc')
ds_ba_years = xr.open_dataset(ruta_gfed4+'GFED4_BA_years.nc')
ds_ba_months = xr.open_dataset(ruta_gfed4+'GFED4_BA_months.nc')
ds_ba_crop = xr.open_dataset(ruta_gfed4+'GFED4_BA_croplands.nc')
BA_a = ds_ba_a.BurnedArea		# Burned Area annal mean
BA_m = ds_ba_m.BurnedArea	# Burned Area monthly means
BA_years = ds_ba_years.BurnedArea		# Burned Area by year
BA_months = ds_ba_months.BurnedArea	# Burned Area by month
BA_crop = ds_ba_crop.BurnedArea 	# Cropland percentage of Burned Area annal mean
BA = BA_years.max(dim='year')

ruta_wfde5='../../DATA/WFDE5/'	# WFDE5 data path
ds_tp_a = xr.open_dataset(ruta_wfde5+'wfde5_tp_a_mean_interp025.nc')
ds_t2m_a = xr.open_dataset(ruta_wfde5+'wfde5_t2m_a_mean_interp025.nc')
ds_tp_m = xr.open_dataset(ruta_wfde5+'wfde5_tp_m_mean_interp025.nc')
ds_t2m_m = xr.open_dataset(ruta_wfde5+'wfde5_t2m_m_mean_interp025.nc')
MAP = ds_tp_a.tp			# Total Precipitation annal mean
MAT = ds_t2m_a.t2m			# 2m Temperature annal mean
Pm = ds_tp_m.tp				# Total Precipitation monthly means
Tm = ds_t2m_m.t2m			# 2m Temperature monthly means
Tmax = Tm.max(dim='month') 	# Temperature of the hotest month
Pmin = Pm.min(dim='month')	# Total precipitation of the driest month

ds_tp_months = xr.open_dataset(ruta_wfde5+'wfde5_tp_months_interp025.nc')
ds_t2m_months = xr.open_dataset(ruta_wfde5+'wfde5_t2m_months_interp025.nc')
P_months = ds_tp_months.tp 										# Month total precipitation
T_months = ds_t2m_months.t2m 									# Month mean temperature
P_years = P_months.groupby('time.year').sum(dim='time')			# Year total precipitation
T_years = T_months.groupby('time.year').mean(dim='time')		# Year mean temperature
Pmin_years = P_months.groupby('time.year').min(dim='time')		# Total precipitation of the driest month by year
Tmax_years = T_months.groupby('time.year').max(dim='time')		# Temperature of the hottest month by year

ds_kg = xr.open_dataset('KG_class.nc')
ds_fs = xr.open_dataset('FS_array.nc')

KG = ds_kg.KG 						# KG classification
FS = ds_fs.FS_array 				# Binary array over 12 months with 1 in fire season months and 0 in non fire season months. Shape: time-12,lat-720,lon-1440

BA_years_m = BA_months.copy()
for i in range(BA_years.year.size):		# Loop over the years
	year = BA_years.year[i].values		# Year value
	for j in range(12):
		BA_years_m.values[j+12*i,:,:] = BA_years.sel(year=year).values[:]

### Some considerations ###
###########################

# We exclude from the analysis locations with >90% of mean annual BA in croplands
P_years = xr.where(BA_crop>90.,np.nan,P_years)
T_years = xr.where(BA_crop>90.,np.nan,T_years)
Pmin_years = xr.where(BA_crop>90.,np.nan,Pmin_years)
Tmax_years = xr.where(BA_crop>90.,np.nan,Tmax_years)
P_months = xr.where(BA_crop>90.,np.nan,P_months)
T_months = xr.where(BA_crop>90.,np.nan,T_months)
BA_years = xr.where(BA_crop>90.,np.nan,BA_years)
BA_months = xr.where(BA_crop>90.,np.nan,BA_months)


# We consider locations as fire-prone if max annual BA>100ha
BA_thres = 100.
# Fire-prone
P_years_p = xr.where(BA_years<BA_thres,np.nan,P_years)
T_years_p = xr.where(BA_years<BA_thres,np.nan,T_years)
Pmin_years_p = xr.where(BA_years<BA_thres,np.nan,Pmin_years)
Tmax_years_p = xr.where(BA_years<BA_thres,np.nan,Tmax_years)
P_months_p = xr.where((BA_years_m<BA_thres),np.nan,P_months)
T_months_p = xr.where((BA_years_m<BA_thres),np.nan,T_months)
# Non fire-prone
P_years_n = xr.where(BA_years>=BA_thres,np.nan,P_years)
T_years_n = xr.where(BA_years>=BA_thres,np.nan,T_years)
Pmin_years_n = xr.where(BA_years>=BA_thres,np.nan,Pmin_years)
Tmax_years_n = xr.where(BA_years>=BA_thres,np.nan,Tmax_years)
P_months_n = xr.where((BA_years_m>=BA_thres),np.nan,P_months)
T_months_n = xr.where((BA_years_m>=BA_thres),np.nan,T_months)

PA_p = P_years_p.mean(dim='year')
TA_p = T_years_p.mean(dim='year')
PMIN_p = Pmin_years_p.mean(dim='year')
TMAX_p = Tmax_years_p.mean(dim='year')
PM_p = P_months_p.groupby('time.month').mean(dim='time')
TM_p = T_months_p.groupby('time.month').mean(dim='time')
PM_p = xr.where(FS==0.,np.nan,PM_p) 	# We exclude months not included in the fire season
TM_p = xr.where(FS==0.,np.nan,TM_p) 	# We exclude months not included in the fire season

PA_n = P_years_n.mean(dim='year')
TA_n = T_years_n.mean(dim='year')
PMIN_n = Pmin_years_n.mean(dim='year')
TMAX_n = Tmax_years_n.mean(dim='year')
PM_n_interannual = P_months_n.groupby('time.month').mean(dim='time')
TM_n_interannual = T_months_n.groupby('time.month').mean(dim='time')
PM_n_season = P_months.groupby('time.month').mean(dim='time')
TM_n_season = T_months.groupby('time.month').mean(dim='time')
PM_n_interannual = xr.where(FS==0.,np.nan,PM_n_interannual) # We exclude months not included in the fire season
TM_n_interannual = xr.where(FS==0.,np.nan,TM_n_interannual) # We exclude months not included in the fire season
PM_n_season = xr.where(FS==1.,np.nan,PM_n_season) 			# We exclude months included in the fire season
TM_n_season = xr.where(FS==1.,np.nan,TM_n_season) 			# We exclude months included in the fire season

### Distribution plots ###
##########################
def distribution(x_p,p_p,x,p_n,KG,var,thres,dx,thres_area,per_p_thres,per_n_thres,baper,decimals):
	if (var=='Pa') | (var=='Pm_season') | (var=='Pm_interannual') | (var=='Pmin'):
		units='mm'
	elif (var=='Ta') | (var=='Tm_season') | (var=='Tm_interannual') | (var=='Tmax'):
		units='$^\circ C$'
	color=['g','y','#ff6bcb','dodgerblue']
	fig = plt.figure(figsize=(8,6))
	ax = fig.add_subplot(111)
	plt.plot(x_p, p_p,color[KG-1])
	plt.plot(x_p, p_n,'k')
	plt.xlabel(var+' ('+units+')')
	plt.ylabel('density')
	ax.axvline(thres,color='k')
	if (var=='Pm_season') | (var=='Pm_interannual') | (var=='Pmin'):
		ax.fill_between(x_p[x_p<=thres], p_p[x_p<=thres], p_n[x_p<=thres], where=(p_p[x_p<=thres] >= p_n[x_p<=thres]), color=color[KG-1], alpha=0.5)
		ax.fill_between(x_p[x_p<=thres], p_p[x_p<=thres], p_n[x_p<=thres], where=(p_p[x_p<=thres] <= p_n[x_p<=thres]), color='grey', alpha=0.5)
		if var=='Pm':
			plt.xlim(-20,1000)
		elif (var=='Pmin') & (KG!=2):
			plt.xlim(-6,150)
	elif (var=='Pa') | (var=='Ta') | (var=='Tm_season') | (var=='Tm_interannual') | (var=='Tmax'):
		ax.fill_between(x_p[x_p>=thres], p_p[x_p>=thres], p_n[x_p>=thres], where=(p_p[x_p>=thres] >= p_n[x_p>=thres]), color=color[KG-1], alpha=0.5)
		ax.fill_between(x_p[x_p>=thres], p_p[x_p>=thres], p_n[x_p>=thres], where=(p_p[x_p>=thres] <= p_n[x_p>=thres]), color='grey', alpha=0.5)
	ax.yaxis.set_major_formatter(FormatStrFormatter('%.1e')) 
	textstr = '\n'.join((
    	'Threshold: {:.{}f} \u00B1 {:.2g} {}'.format(thres,decimals,dx,units),
    	'Area difference: {:.2f}'.format(thres_area),
    	'Fire-prone classified: {:.1f}%'.format(float(per_p_thres)),
       	'Non fire-prone classified: {:.1f}%'.format(float(per_n_thres)),
    	'Corresponding BA: {:.1f}%'.format(float(baper))))
	props = dict(boxstyle='round', facecolor='w', alpha=0.75)
	if (var=='Pm_season') | (var=='Pm_interannual') | (var=='Pmin') | (var=='Pa'):
		pos_x = 0.5
	else:
		pos_x = 0.05
	ax.text(pos_x, 0.95, textstr, transform=ax.transAxes, fontsize=12,
        verticalalignment='top', horizontalalignment='left',bbox=props)
	plt.grid(True)
	ratio=0.8
	ax.set_aspect(1.0/ax.get_data_ratio()*ratio)
	fig.savefig('./Distribution_plots/'+var+'_KG'+str(KG)+'_dist.png',dpi=150)


### Distributions ###
#####################

# Dataframe for saving the selected thresholds
df = pd.DataFrame({'FC': [1, 2, 3, 4],
					'Pa': [np.nan, np.nan, np.nan, np.nan],
					'Ta': [np.nan, np.nan, np.nan, np.nan],
					'Pmin': [np.nan, np.nan, np.nan, np.nan],
					'Tmax': [np.nan, np.nan, np.nan, np.nan],
					'Pms': [np.nan, np.nan, np.nan, np.nan],
					'Tms': [np.nan, np.nan, np.nan, np.nan],
					'Pmi': [np.nan, np.nan, np.nan, np.nan],
					'Tmi': [np.nan, np.nan, np.nan, np.nan]})
# Dataframe for saving the selected thresholds uncertainty
df_dx = pd.DataFrame({'FC': [1, 2, 3, 4],
					'Pa': [np.nan, np.nan, np.nan, np.nan],
					'Ta': [np.nan, np.nan, np.nan, np.nan],
					'Pmin': [np.nan, np.nan, np.nan, np.nan],
					'Tmax': [np.nan, np.nan, np.nan, np.nan],
					'Pms': [np.nan, np.nan, np.nan, np.nan],
					'Tms': [np.nan, np.nan, np.nan, np.nan],
					'Pmi': [np.nan, np.nan, np.nan, np.nan],
					'Tmi': [np.nan, np.nan, np.nan, np.nan]})


# Loop over Koppen-Geiger classes
for kg_class in range(1,5):

	Pa_p_a = PA_p.copy()
	Pa_n_a = PA_n.copy()
	Pa_p_a = xr.where((KG!=kg_class),np.nan,Pa_p_a)
	Pa_n_a = xr.where((KG!=kg_class),np.nan,Pa_n_a)
	Pa_p = Pa_p_a.values[np.logical_not(np.isnan(Pa_p_a.values))]
	Pa_n = Pa_n_a.values[np.logical_not(np.isnan(Pa_n_a.values))]

	Ta_p_a = TA_p.copy()
	Ta_n_a = TA_n.copy()
	Ta_p_a = xr.where((KG!=kg_class),np.nan,Ta_p_a)
	Ta_n_a = xr.where((KG!=kg_class),np.nan,Ta_n_a)
	Ta_p = Ta_p_a.values[np.logical_not(np.isnan(Ta_p_a.values))]
	Ta_n = Ta_n_a.values[np.logical_not(np.isnan(Ta_n_a.values))]

	Tmax_p_a = TMAX_p.copy()
	Tmax_n_a = TMAX_n.copy()
	Tmax_p_a = xr.where((KG!=kg_class),np.nan,Tmax_p_a)
	Tmax_n_a = xr.where((KG!=kg_class),np.nan,Tmax_n_a)
	Tmax_p = Tmax_p_a.values[np.logical_not(np.isnan(Tmax_p_a.values))]
	Tmax_n = Tmax_n_a.values[np.logical_not(np.isnan(Tmax_n_a.values))]

	Pmin_p_a = PMIN_p.copy()
	Pmin_n_a = PMIN_n.copy()
	Pmin_p_a = xr.where((KG!=kg_class),np.nan,Pmin_p_a)
	Pmin_n_a = xr.where((KG!=kg_class),np.nan,Pmin_n_a)
	Pmin_p = Pmin_p_a.values[np.logical_not(np.isnan(Pmin_p_a.values))]
	Pmin_n = Pmin_n_a.values[np.logical_not(np.isnan(Pmin_n_a.values))]

	Pm_p_a = PM_p.copy()
	Pm_p_a = xr.where((KG!=kg_class),np.nan,Pm_p_a)
	Pm_p = Pm_p_a.values[np.logical_not(np.isnan(Pm_p_a.values))]
	Pm_n_a = PM_n_season.copy()
	Pm_n_a = xr.where((KG!=kg_class),np.nan,Pm_n_a)
	Pm_n_s = Pm_n_a.values[np.logical_not(np.isnan(Pm_n_a.values))]
	Pm_n_a = PM_n_interannual.copy()
	Pm_n_a = xr.where((KG!=kg_class),np.nan,Pm_n_a)
	Pm_n_i = Pm_n_a.values[np.logical_not(np.isnan(Pm_n_a.values))]

	Tm_p_a = TM_p.copy()
	Tm_p_a = xr.where((KG!=kg_class),np.nan,Tm_p_a)
	Tm_p = Tm_p_a.values[np.logical_not(np.isnan(Tm_p_a.values))]
	Tm_n_a = TM_n_season.copy()
	Tm_n_a = xr.where((KG!=kg_class),np.nan,Tm_n_a)
	Tm_n_s = Tm_n_a.values[np.logical_not(np.isnan(Tm_n_a.values))]
	Tm_n_a = TM_n_interannual.copy()
	Tm_n_a = xr.where((KG!=kg_class),np.nan,Tm_n_a)
	Tm_n_i = Tm_n_a.values[np.logical_not(np.isnan(Tm_n_a.values))]

	BA_years_kg = xr.where(KG!=kg_class,0.,BA_years)
	BA_months_kg = xr.where(KG!=kg_class,0.,BA_months)


	### Pa ###
	range_min = np.min([Pa_p.min(),Pa_n.min()])
	range_max = np.max([Pa_p.max(),Pa_n.max()])
	dx = 15.
	n = int((range_max-range_min)/dx)		# n is the number of bins (variable discretization)
	# We obtain the probability density function for the fire prone data (p_p) and for the non fire-prone data (p_n)
	p_p, x_p = np.histogram(Pa_p, bins=n, range=(range_min,range_max), density=True)
	p_n, x_n = np.histogram(Pa_n, bins=n, range=(range_min,range_max), density=True)
	x = x_p[:-1]
	x_p = x_p[:-1] + (x_p[1] - x_p[0])/2   # convert bin edges to centers
	dx_map = x_p[1]-x_p[0]
	# We obtain the area of the density function enclosed by different thresholds for the fire prone data and for the non fire-prone data
	# The greater the difference of these two areas, the greater the threshold performance is
	dif_area = np.zeros(p_p.size)
	for i in range(0,p_p.size):
		dif_area[i] = (p_p[x[:]>=x[i]].sum()*dx_map)-(p_n[x[:]>=x[i]].sum()*dx_map)
	# The variable threshold is the one that maximizes the area difference. More fire-prone points, less non fire-prone points.
	map_thres=x[np.argmax(dif_area)]
	map_thres_area = np.max(dif_area)
	# We compute the percentage of fire-prone points that meet the threshold (map_per_p_thres)
	map_per_p_thres = (p_p[x[:]>=map_thres].sum()*dx_map)*100.
	# We compute the percentage of non fire-prone points that meet the threshold (map_per_n_thres)
	map_per_n_thres = (p_n[x[:]>=map_thres].sum()*dx_map)*100.
	# We compute the percentage of burned area corresponding to the fire-prone points that meet the threshold (map_baper)
	map_baper = BA_years_kg.where(P_years>map_thres).sum()/BA_years_kg.sum()*100.
	# Probability density function plot
	dx_map = dx_map*0.5
	decimals_map = '{:.2g}'.format(dx_map)[::-1].find('.')
	if decimals_map==-1:
		decimals_map=0
		map_thres = round(map_thres)
		if len('{:.2g}'.format(dx_map))==1:
			decimals_map=1
	else:
		map_thres = round(map_thres,decimals_map)
	distribution(x_p,p_p,x,p_n,kg_class,'Pa',map_thres,dx_map,map_thres_area,map_per_p_thres,map_per_n_thres,map_baper,decimals_map)
	# We do the same for each variable

	### Ta ###
	range_min = np.min([Ta_p.min(),Ta_n.min()])
	range_max = np.max([Ta_p.max(),Ta_n.max()])
	dx = 0.2
	n = int((range_max-range_min)/dx)		# n is the number of bins (variable discretization)
	p_p, x_p = np.histogram(Ta_p, bins=n, range=(range_min,range_max), density=True)
	p_n, x_n = np.histogram(Ta_n, bins=n, range=(range_min,range_max), density=True)
	x = x_p[:-1]
	x_p = x_p[:-1] + (x_p[1] - x_p[0])/2   # convert bin edges to centers
	dx_mat = x_p[1]-x_p[0]
	dif_area = np.zeros(p_p.size)
	for i in range(0,p_p.size):
		dif_area[i] = (p_p[x[:]>=x[i]].sum()*dx_mat)-(p_n[x[:]>=x[i]].sum()*dx_mat)
	mat_thres=x[np.argmax(dif_area)]
	mat_thres_area = np.max(dif_area)
	mat_per_p_thres = (p_p[x[:]>=mat_thres].sum()*dx_mat)*100.
	mat_per_n_thres = (p_n[x[:]>=mat_thres].sum()*dx_mat)*100.
	mat_baper = BA_years_kg.where(T_years>mat_thres).sum()/BA_years_kg.sum()*100.
	dx_mat = dx_mat*0.5
	decimals_mat = '{:.2g}'.format(dx_mat)[::-1].find('.')
	if decimals_mat==-1:
		decimals_mat=0
		mat_thres = round(mat_thres)
		if len('{:.2g}'.format(dx_mat))==1:
			decimals_mat=1
	else:
		mat_thres = round(mat_thres,decimals_mat)
	distribution(x_p,p_p,x,p_n,kg_class,'Ta',mat_thres,dx_mat,mat_thres_area,mat_per_p_thres,mat_per_n_thres,mat_baper,decimals_mat)

	### Pmin ###
	range_min = np.min([Pmin_p.min(),Pmin_n.min()])
	range_max = np.max([Pmin_p.max(),Pmin_n.max()])
	dx = 1.
	n = int((range_max-range_min)/dx)		# n is the number of bins (variable discretization)
	p_p, x_p = np.histogram(Pmin_p, bins=n, range=(range_min,range_max), density=True)
	p_n, x_n = np.histogram(Pmin_n, bins=n, range=(range_min,range_max), density=True)
	x = x_p[1:]
	x_p = x_p[:-1] + (x_p[1] - x_p[0])/2   # convert bin edges to centers
	dx_pmin = x_p[1]-x_p[0]
	dif_area = np.zeros(p_p.size)
	for i in range(0,p_p.size):
		dif_area[i] = (p_p[x[:]<=x[i]].sum()*dx_pmin)-(p_n[x[:]<=x[i]].sum()*dx_pmin)
	pmin_thres=x[np.argmax(dif_area)]
	pmin_thres_area = np.max(dif_area)
	pmin_per_p_thres = (p_p[x[:]<=pmin_thres].sum()*dx_pmin)*100.
	pmin_per_n_thres = (p_n[x[:]<=pmin_thres].sum()*dx_pmin)*100.
	pmin_baper = BA_years_kg.where(Pmin_years<pmin_thres).sum()/BA_years_kg.sum()*100.
	dx_pmin = dx_pmin*0.5
	decimals_pmin = '{:.2g}'.format(dx_pmin)[::-1].find('.')
	if decimals_pmin==-1:
		decimals_pmin=0
		pmin_thres = round(pmin_thres)
		if len('{:.2g}'.format(dx_pmin))==1:
			decimals_pmin=1
	else:
		pmin_thres = round(pmin_thres,decimals_pmin)
	distribution(x_p,p_p,x,p_n,kg_class,'Pmin',pmin_thres,dx_pmin,pmin_thres_area,pmin_per_p_thres,pmin_per_n_thres,pmin_baper,decimals_pmin)

	### Tmax ###
	range_min = np.min([Tmax_p.min(),Tmax_n.min()])
	range_max = np.max([Tmax_p.max(),Tmax_n.max()])
	dx = 0.3
	n = int((range_max-range_min)/dx)		# n is the number of bins (variable discretization)
	p_p, x_p = np.histogram(Tmax_p, bins=n, range=(range_min,range_max), density=True)
	p_n, x_n = np.histogram(Tmax_n, bins=n, range=(range_min,range_max), density=True)
	x = x_p[:-1]
	x_p = x_p[:-1] + (x_p[1] - x_p[0])/2   # convert bin edges to centers
	dx_tmax = x_p[1]-x_p[0]
	dif_area = np.zeros(p_p.size)
	for i in range(0,p_p.size):
		dif_area[i] = (p_p[x[:]>=x[i]].sum()*dx_tmax)-(p_n[x[:]>=x[i]].sum()*dx_tmax)
	tmax_thres=x[np.argmax(dif_area)]
	tmax_thres_area = np.max(dif_area)
	tmax_per_p_thres = (p_p[x[:]>=tmax_thres].sum()*dx_tmax)*100.
	tmax_per_n_thres = (p_n[x[:]>=tmax_thres].sum()*dx_tmax)*100.
	tmax_baper = BA_years_kg.where(Tmax_years>tmax_thres).sum()/BA_years_kg.sum()*100.
	dx_tmax = dx_tmax*0.5
	decimals_tmax = '{:.2g}'.format(dx_tmax)[::-1].find('.')
	if decimals_tmax==-1:
		decimals_tmax=0
		tmax_thres = round(tmax_thres)
		if len('{:.2g}'.format(dx_tmax))==1:
			decimals_tmax=1
	else:
		tmax_thres = round(tmax_thres,decimals_tmax)
	distribution(x_p,p_p,x,p_n,kg_class,'Tmax',tmax_thres,dx_tmax,tmax_thres_area,tmax_per_p_thres,tmax_per_n_thres,tmax_baper,decimals_tmax)

	### Pm ###
	range_min = np.min([Pm_p.min(),Pm_n_s.min()])
	range_max = np.max([Pm_p.max(),Pm_n_s.max()])
	dx = 3.
	n = int((range_max-range_min)/dx)		# n is the number of bins (variable discretization)
	p_p, x_p = np.histogram(Pm_p, bins=n, range=(range_min,range_max), density=True)
	p_n, x_n = np.histogram(Pm_n_s, bins=n, range=(range_min,range_max), density=True)
	x = x_p[1:]
	x_p = x_p[:-1] + (x_p[1] - x_p[0])/2   # convert bin edges to centers
	dx_pms = x_p[1]-x_p[0]
	dif_area = np.zeros(p_p.size)
	for i in range(0,p_p.size):
		dif_area[i] = (p_p[x[:]<=x[i]].sum()*dx_pms)-(p_n[x[:]<=x[i]].sum()*dx_pms)
	pms_thres=x[np.argmax(dif_area)]
	pms_thres_area = np.max(dif_area)
	pms_per_p_thres = (p_p[x[:]<=pms_thres].sum()*dx_pms)*100.
	pms_per_n_thres = (p_n[x[:]<=pms_thres].sum()*dx_pms)*100.
	pms_baper = BA_months_kg.where(P_months<pms_thres).sum()/BA_months_kg.sum()*100.
	dx_pms = dx_pms*0.5
	decimals_pms = '{:.2g}'.format(dx_pms)[::-1].find('.')
	if decimals_pms==-1:
		decimals_pms=0
		pms_thres = round(pms_thres)
		if len('{:.2g}'.format(dx_pms))==1:
			decimals_pms=1
	else:
		pms_thres = round(pms_thres,decimals_pms)
	distribution(x_p,p_p,x,p_n,kg_class,'Pm_season',pms_thres,dx_pms,pms_thres_area,pms_per_p_thres,pms_per_n_thres,pms_baper,decimals_pms)

	### Tm ###
	range_min = np.min([Tm_p.min(),Tm_n_s.min()])
	range_max = np.max([Tm_p.max(),Tm_n_s.max()])
	dx = 0.5
	n = int((range_max-range_min)/dx)		# n is the number of bins (variable discretization)
	p_p, x_p = np.histogram(Tm_p, bins=n, range=(range_min,range_max), density=True)
	p_n, x_n = np.histogram(Tm_n_s, bins=n, range=(range_min,range_max), density=True)
	x = x_p[:-1]
	x_p = x_p[:-1] + (x_p[1] - x_p[0])/2   # convert bin edges to centers
	dx_tms = x_p[1]-x_p[0]
	dif_area = np.zeros(p_p.size)
	for i in range(0,p_p.size):
		dif_area[i] = (p_p[x[:]>=x[i]].sum()*dx_tms)-(p_n[x[:]>=x[i]].sum()*dx_tms)
	tms_thres=x[np.argmax(dif_area)]
	tms_thres_area = np.max(dif_area)
	tms_per_p_thres = (p_p[x[:]>=tms_thres].sum()*dx_tms)*100.
	tms_per_n_thres = (p_n[x[:]>=tms_thres].sum()*dx_tms)*100.
	tms_baper = BA_months_kg.where(T_months>tms_thres).sum()/BA_months_kg.sum()*100.
	dx_tms = dx_tms*0.5
	decimals_tms = '{:.2g}'.format(dx_tms)[::-1].find('.')
	if decimals_tms==-1:
		decimals_tms=0
		tms_thres = round(tms_thres)
		if len('{:.2g}'.format(dx_tms))==1:
			decimals_tms=1
	else:
		tms_thres = round(tms_thres,decimals_tms)
	distribution(x_p,p_p,x,p_n,kg_class,'Tm_season',tms_thres,dx_tms,tms_thres_area,tms_per_p_thres,tms_per_n_thres,tms_baper,decimals_tms)

	### Pm ###
	range_min = np.min([Pm_p.min(),Pm_n_i.min()])
	range_max = np.max([Pm_p.max(),Pm_n_i.max()])
	dx = 2.
	n = int((range_max-range_min)/dx)		# n is the number of bins (variable discretization)
	p_p, x_p = np.histogram(Pm_p, bins=n, range=(range_min,range_max), density=True)
	p_n, x_n = np.histogram(Pm_n_i, bins=n, range=(range_min,range_max), density=True)
	x = x_p[1:]
	x_p = x_p[:-1] + (x_p[1] - x_p[0])/2   # convert bin edges to centers
	dx_pmi = x_p[1]-x_p[0]
	dif_area = np.zeros(p_p.size)
	for i in range(0,p_p.size):
		dif_area[i] = (p_p[x[:]<=x[i]].sum()*dx_pmi)-(p_n[x[:]<=x[i]].sum()*dx_pmi)
	pmi_thres=x[np.argmax(dif_area)]
	pmi_thres_area = np.max(dif_area)
	pmi_per_p_thres = (p_p[x[:]<=pmi_thres].sum()*dx_pmi)*100.
	pmi_per_n_thres = (p_n[x[:]<=pmi_thres].sum()*dx_pmi)*100.
	pmi_baper = BA_months_kg.where(P_months<pmi_thres).sum()/BA_months_kg.sum()*100.
	dx_pmi = dx_pmi*0.5
	decimals_pmi = '{:.2g}'.format(dx_pmi)[::-1].find('.')
	if decimals_pmi==-1:
		decimals_pmi=0
		pmi_thres = round(pmi_thres)
		if len('{:.2g}'.format(dx_pmi))==1:
			decimals_pmi=1
	else:
		pmi_thres = round(pmi_thres,decimals_pmi)
	distribution(x_p,p_p,x,p_n,kg_class,'Pm_interannual',pmi_thres,dx_pmi,pmi_thres_area,pmi_per_p_thres,pmi_per_n_thres,pmi_baper,decimals_pmi)

	### Tm ###
	range_min = np.min([Tm_p.min(),Tm_n_i.min()])
	range_max = np.max([Tm_p.max(),Tm_n_i.max()])
	dx = 0.5
	n = int((range_max-range_min)/dx)		# n is the number of bins (variable discretization)
	p_p, x_p = np.histogram(Tm_p, bins=n, range=(range_min,range_max), density=True)
	p_n, x_n = np.histogram(Tm_n_i, bins=n, range=(range_min,range_max), density=True)
	x = x_p[:-1]
	x_p = x_p[:-1] + (x_p[1] - x_p[0])/2   # convert bin edges to centers
	dx_tmi = x_p[1]-x_p[0]
	dif_area = np.zeros(p_p.size)
	for i in range(0,p_p.size):
		dif_area[i] = (p_p[x[:]>=x[i]].sum()*dx_tmi)-(p_n[x[:]>=x[i]].sum()*dx_tmi)
	tmi_thres=x[np.argmax(dif_area)]
	tmi_thres_area = np.max(dif_area)
	tmi_per_p_thres = (p_p[x[:]>=tmi_thres].sum()*dx_tmi)*100.
	tmi_per_n_thres = (p_n[x[:]>=tmi_thres].sum()*dx_tmi)*100.
	tmi_baper = BA_months_kg.where(T_months>tmi_thres).sum()/BA_months_kg.sum()*100.
	dx_tmi = dx_tmi*0.5
	decimals_tmi = '{:.2g}'.format(dx_tmi)[::-1].find('.')
	if decimals_tmi==-1:
		decimals_tmi=0
		tmi_thres = round(tmi_thres)
		if len('{:.2g}'.format(dx_tmi))==1:
			decimals_tmi=1
	else:
		tmi_thres = round(tmi_thres,decimals_tmi)
	distribution(x_p,p_p,x,p_n,kg_class,'Tm_interannual',tmi_thres,dx_tmi,tmi_thres_area,tmi_per_p_thres,tmi_per_n_thres,tmi_baper,decimals_tmi)

	# We print the data
	print('\n')
	print('KG{}, Pa threshold: {:.{}f} \u00B1 {:.{}f}'.format(kg_class,map_thres,decimals_map,dx_map,decimals_map))
	print('KG{}, Pa area difference: {:.3f}'.format(kg_class,map_thres_area))
	print('KG{}, Pa classified fire-prone points: {:.1f}%'.format(kg_class,float(map_per_p_thres)))
	print('KG{}, Pa classified non-fire-prone points: {:.1f}%'.format(kg_class,float(map_per_n_thres)))
	print('KG{}, Pa BA of classified points: {:.1f}%'.format(kg_class,float(map_baper)))
	print('KG{}, Ta threshold: {:.{}f} \u00B1 {:.{}f}'.format(kg_class,mat_thres,decimals_mat,dx_mat,decimals_mat))
	print('KG{}, Ta area difference: {:.3f}'.format(kg_class,mat_thres_area))
	print('KG{}, Ta classified fire-prone points: {:.1f}%'.format(kg_class,float(mat_per_p_thres)))
	print('KG{}, Ta classified non-fire-prone points: {:.1f}%'.format(kg_class,float(mat_per_n_thres)))
	print('KG{}, Ta BA of classified points: {:.1f}%'.format(kg_class,float(mat_baper)))
	print('KG{}, Pmin threshold: {:.{}f} \u00B1 {:.{}f}'.format(kg_class,pmin_thres,decimals_pmin,dx_pmin,decimals_pmin))
	print('KG{}, Pmin area difference: {:.3f}'.format(kg_class,pmin_thres_area))
	print('KG{}, Pmin classified fire-prone points: {:.1f}%'.format(kg_class,float(pmin_per_p_thres)))
	print('KG{}, Pmin classified non-fire-prone points: {:.1f}%'.format(kg_class,float(pmin_per_n_thres)))
	print('KG{}, Pmin BA of classified points: {:.1f}%'.format(kg_class,float(pmin_baper)))
	print('KG{}, Tmax threshold: {:.{}f} \u00B1 {:.{}f}'.format(kg_class,tmax_thres,decimals_tmax,dx_tmax,decimals_tmax))
	print('KG{}, Tmax area difference: {:.3f}'.format(kg_class,tmax_thres_area))
	print('KG{}, Tmax classified fire-prone points: {:.1f}%'.format(kg_class,float(tmax_per_p_thres)))
	print('KG{}, Tmax classified non-fire-prone points: {:.1f}%'.format(kg_class,float(tmax_per_n_thres)))
	print('KG{}, Tmax BA of classified points: {:.1f}%'.format(kg_class,float(tmax_baper)))
	print('KG{}, Pm_season threshold: {:.{}f} \u00B1 {:.{}f}'.format(kg_class,pms_thres,decimals_pms,dx_pms,decimals_pms))
	print('KG{}, Pm_season area difference: {:.3f}'.format(kg_class,pms_thres_area))
	print('KG{}, Pm_season classified fire-prone points: {:.1f}%'.format(kg_class,float(pms_per_p_thres)))
	print('KG{}, Pm_season classified non-fire-prone points: {:.1f}%'.format(kg_class,float(pms_per_n_thres)))
	print('KG{}, Pm_season BA of classified points: {:.1f}%'.format(kg_class,float(pms_baper)))
	print('KG{}, Tm_season threshold: {:.{}f} \u00B1 {:.{}f}'.format(kg_class,tms_thres,decimals_tms,dx_tms,decimals_tms))
	print('KG{}, Tm_season area difference: {:.3f}'.format(kg_class,tms_thres_area))
	print('KG{}, Tm_season classified fire-prone points: {:.1f}%'.format(kg_class,float(tms_per_p_thres)))
	print('KG{}, Tm_season classified non-fire-prone points: {:.1f}%'.format(kg_class,float(tms_per_n_thres)))
	print('KG{}, Tm_season BA of classified points: {:.1f}%'.format(kg_class,float(tms_baper)))
	print('KG{}, Pm_interannual threshold: {:.{}f} \u00B1 {:.{}f}'.format(kg_class,pmi_thres,decimals_pmi,dx_pmi,decimals_pmi))
	print('KG{}, Pm_interannual area difference: {:.3f}'.format(kg_class,pmi_thres_area))
	print('KG{}, Pm_interannual classified fire-prone points: {:.1f}%'.format(kg_class,float(pmi_per_p_thres)))
	print('KG{}, Pm_interannual classified non-fire-prone points: {:.1f}%'.format(kg_class,float(pmi_per_n_thres)))
	print('KG{}, Pm_interannual BA of classified points: {:.1f}%'.format(kg_class,float(pmi_baper)))
	print('KG{}, Tm_interannual threshold: {:.{}f} \u00B1 {:.{}f}'.format(kg_class,tmi_thres,decimals_tmi,dx_tmi,decimals_tmi))
	print('KG{}, Tm_interannual area difference: {:.3f}'.format(kg_class,tmi_thres_area))
	print('KG{}, Tm_interannual classified fire-prone points: {:.1f}%'.format(kg_class,float(tmi_per_p_thres)))
	print('KG{}, Tm_interannual classified non-fire-prone points: {:.1f}%'.format(kg_class,float(tmi_per_n_thres)))
	print('KG{}, Tm_interannual BA of classified points: {:.1f}%'.format(kg_class,float(tmi_baper)))

	# Between the 6 variables, we select the those that have better results in terms of the density function area difference
	var_number = np.array([0,1,2,3,4,5,6,7])
	thres = np.array([map_thres,mat_thres,pmin_thres,tmax_thres,pms_thres,tms_thres,pmi_thres,tmi_thres])
	dx_thres = np.array([dx_map,dx_mat,dx_pmin,dx_tmax,dx_pms,dx_tms,dx_pmi,dx_tmi])
	thres_area = np.array([map_thres_area,mat_thres_area,pmin_thres_area,tmax_thres_area,pms_thres_area,tms_thres_area,pmi_thres_area,tmi_thres_area])
	p_thres_area = np.array([map_per_p_thres,mat_per_p_thres,pmin_per_p_thres,tmax_per_p_thres,pms_per_p_thres,tms_per_p_thres,pmi_per_p_thres,tmi_per_p_thres])

	# We save the threshold in the data frame
	for i in var_number:
		df.iloc[kg_class-1,var_number+1]=thres[var_number]
		df_dx.iloc[kg_class-1,var_number+1]=dx_thres[var_number] 
	if (kg_class==2):
		df.iloc[kg_class-1,5]=round(np.percentile(Pm_p,90),1)	
		df_dx.iloc[kg_class-1,5]=0.
	if (kg_class==4):
		df.iloc[kg_class-1,7]=round(np.percentile(Pm_p,75),1)	
		df_dx.iloc[kg_class-1,5]=0.
		
# We use the Arid Pa threshold for those climates that has not a Pa threshold
for i in [0,2,3]:
	df.loc[int(i),'Pa']=df.loc[1,'Pa']
	df_dx.loc[int(i),'Pa']=df_dx.loc[1,'Pa']

# Saving the threshold data frame as a csv file
df.to_csv('./Thresholds.csv')
df_dx['Pa'] = df_dx['Pa'].astype('float32').map('{:.2g}'.format)
df_dx['Ta'] = df_dx['Ta'].astype('float32').map('{:.2g}'.format)
df_dx['Pmin'] = df_dx['Pmin'].astype('float32').map('{:.2g}'.format)
df_dx['Tmax'] = df_dx['Tmax'].astype('float32').map('{:.2g}'.format)
df_dx['Pms'] = df_dx['Pms'].astype('float32').map('{:.2g}'.format)
df_dx['Tms'] = df_dx['Tms'].astype('float32').map('{:.2g}'.format)
df_dx['Pmi'] = df_dx['Pmi'].astype('float32').map('{:.2g}'.format)
df_dx['Tmi'] = df_dx['Tmi'].astype('float32').map('{:.2g}'.format)
df_dx.to_csv('./Thresholds_dx.csv')